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The Ha Variations of the Luminous Blue Variable P Cygni: 
Discrete Absorption Components and the Short S Doradus 

Phase 

N. D. Richardson^ N. D. Morrison^, D. R. Gies\ N. Markova^ E. N. Hesselbach^, 

and J. R. Percy^ 

ABSTRACT 

P Cygni is a prototype of the Luminous Blue Variables (or S Doradus 
variables), and the star displays photometric and emission line variability 
on a timescale of years (known as the "short S Doradus phase" varia- 
tions). Here we present new high resolution Ha spectroscopy of P Cyg 
that we combine with earlier spectra and concurrent V^-band photome- 
try to document the emission and continuum flux variations over a 24 y 
time span. We show that the emission and continuum fluxes vary in con- 
cert on timescales of 1.6 y and longer, but differ on shorter timescales. 
The Ha profile shape also varies on the photometric timescales, and we 
describe the observed co-variations of the emission peak and absorption 
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trough properties. We argue that the episodes of photometric and emis- 
sion brightening are caused by increases in the size of the emission region 
that are related to variations in wind mass loss rate and outflow speed. 
We find evidence of blueward accelerating, Discrete Absorption Compo- 
nents (DACs) in the absorption trough of the Ha profile, and these fea- 
tures have slower accelerations and longer durations than those observed 
in other lines. The DAC strengths also appear to vary on the photomet- 
ric timescales, and we suggest that the propagation of the DAC-related 
wind structures is closely related to changes in the overall wind mass loss 
rate and velocity. 

Subject headings: stars: variables — stars: early-type — stars: individual 
(P Cyg, HD 193237) — stars: winds, outflows — stars: circumstellar 
matter — stars: mass loss 



Introduction 



Luminous Blue Variables (LBVs or S Doradus variables) are evolved, massive 
stars. LBVs are characterized by large mass loss rates and variability on multiple 
timescales. The two "prototypical" Galactic LBVs are r] Carinae and P Cygni, and 
they prob ably represent different extr emes of mass loss rate within the scheme of LBV 
evolution ( llsraelian fc de Grootlll999l ) . One of the defining criteria of the LBVs is the 
observation of a large scale eruption, when the star brightens by several magnitudes. 
The quiescent times between these eruptions may last centuries. In addition to such 
rare, giant eruptions, these stars also display lesser photometric an d spectroscopic 
variations on other timescales (e.g. Humphreys & Davidson 1994). Ivan Genderen 



(l200l[ ) defines the S Doradus (SD-) phase to be the moderate, long-term, brightening 
and fading phases. There are two types of these phases, short and long, with similar 
characteristics. The short SD-phase is typically on the order of years (< 10 years), 
while the longer SD-phase is on a timescale of decades. These phases are thought 
to originate from changes in the star's photosphere, and both may have the same 
physical driving mechanism. These long-term v ariations are observ ed to differ from 
cycle to cycle, both in duration and amplitude ( ISterken et al.lll997l ). 



P Cygni (HD 193237, HR 7763, Nova Cyg 1600) remains one of the most fasci- 
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nating objects in the sky. It was discovered during its first recorded great eruption 
in 1600 by Willem Janszoon Blaeu, a Dutch chart-maker and mathematician. Dur- 
ing this eruption, the star brightened to about 3rd magnitude for about six years, 
and then it faded from visibihty by 1626. It rose again in 1654 to about the same 
maximal brightness, where it remained for five years. The star faded after this, and 
although its long-term variability is poorly docur nented, the star has been slowly 
brightening to its current magnitude of about 4.8 (jlsraelian fc de Grootlll999l ). The 
slow brightening may refiect evolutionary changes (de Groot & Lamers 1992; Lamers 
& de Groot 1992; Langer et al. 1994). 



Lo ng-term photometri c mon itoring of P Cygn i began in the 1980s when lPercy fc Welch 



(Il983l ) , iPercy et al.l (119881 ) , and Ide GrootI (Il990l ) embarked on extended observing 
campaigns. Their observations showed that the variations often occur on three char- 
acteristic timescales: a short ~ 17 day variation similar to the a Cygni type variations 
observed in hot supergiants, a ~ 100 day "quasi-period" similar to that observed in 
other LBVs, and a long-term cycle (years) attributed to a short-SD phase (de Groot 
et al. 2001; Percy et al. 2001). 



According to llsraelian fc de GrootI (119991 ). comprehensive spectral monitoring 
of P Cygni was started by Luud (1967) and Markova (1993), among others. The first 
long-term spectroscopic monitoring campaign of P Cygni was presented in seminal 
papers by Markova et al. (2001a,b). They found evidence of co- variability of the Ha 
emission line strength and Johnson UBV photometry, indicating a short-SD phase 
with a quasi-period of ~7 years, although their observations did not fully cover two 
cycles. The variations were attributed to inversely correlated changes in effective 
temperature and r adius, maintaini n g a ii early constant luminosity. A similar cycle 
time was found by Ide Groot et al.l (l200ll ). who reported on photometric variations 
which were consistent with a timescale of 5.5 to 8.5 years. 

In addition to the large scale variations in emission strength, Markova (2000) 
found that there are at least four other kinds of line profile variability in the spec- 
trum of P Cygni. The most striking of these is the long documented appearance 
of blueward-migrating, absorption sub-features that are called Discrete Absorption 
Components (DACs: Israelian & de Groot 1999; Markova 2000). These are generally 
observed in low and intermediate excitation state lines in the optical (Markova 2000) 
and UV spectrum (Israelian et al. 1996). They are frequently detected in the upper 
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sequence of the H Balmer lines (principal quantum number 9 < n < 15; Markova 
2000), but to our knowledge, DACs have not been reported before now for the absorp- 
tion component of Ha. DACs are often (but not always) narrow (FWHM 10 — 15 
km s~^) and may be unresolved in low dispersion spectra. The DACs tend to appear 
over a radial velocity range of —90 to —200 km s~^ with an acceleration of —0.1 to 
-0.6 km s"^ d-\ A recurrence timescale of ~ 200 d is sometimes observed (Kolka 
1983; Markova 1986a; Israelian et al. 1996; Kolka 1998). These accelerations are 
much slower and the timescales are much longer than those associated with DACs in 
the winds of 0-stars (Kaper et al. 1999). The DACs in the spectrum of P Cygni may 
form in outward moving and dense shells (Kolka 1983; Lamers et al. 1985; Markova 
1986a; Israelian et al. 1996), in spiral-shaped co-rotating interaction regions (CIRs; 
Cranmer & Owocki 1996; Markova 2000), or in dense clumps in the wind (Lepine & 
Moffat 2008). 

In this paper, we present new hi gh resolution Ha spectr oscopy, which we com- 



bined with previous measurements by iMarkova et al.l (l2001al ) to explore the charac- 



teristics of P Cygni's short SD-phase. We also compare this with archival Johnson V 
photometry and new observations obtained by AAVSO observers. Section 2 describes 
our observations. In Section 3, we present our analysis of long-term variations of the 
continuum and the Ha equivalent width. We describe the Ha profile morphology 
changes and DAC propagations in Section 4. Our discussion and conclusions are 
presented in Section 5. 



2. Observations 



We obtained 126 new spectroscopic observation s of P Cygni using t he Ritter 
Observatory 1 m telescope and echelle spectrograph (IMorrison et al.l 119971 ) between 
1999 June 7 and 2007 October 30. These high resolving power {R = 26, 000) spectra 
were reduced by standard techniques with I RAf|^. Observati o ns co llected prior to 
2007 were taken using the setup described in iMorrison et al.l ( ll997l ). These obser- 



^IRAF is distributed by the National Optical Astronomy Observatory, which is operated by the 
Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the 
National Science Foundation. 
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vations record a 70 A range in the order centered on Ha, and they typically have 
a signal-to-noise ratio between 50 and 100 per resolution element in the continuum. 
Observations made during the calendar year 2007 were taken with the same spec- 
trograph, except the camera was a Spectral Instruments 600 Series camera, with a 
front-illuminated Imager Labs IL-C2004 4100x4096 pixel sensor (15x15 micron pix- 
els). To maintain consistency with the older observations, the camera was operated 
with 2x2 pixel binning. The newer observations recorded a larger portion of the 
order centered on Ha and typically reached a signal to noise ratio between 50 and 
100 per resolution element in the continuum. We trimmed these spectra so that 
the wavelength range was the same as in the older data. The spectra taken after 
2002 September have poor wavelength calibration due to problems with the Th-Ar 
lamp. In order to use these spectra for kinematical measurements, the telluric H2O 
lines in the vicinity of Ha were fitted to improve the solution. This worked for most 
cases, but the errors associated with the telluric re-calibration are roughly ±3 km 
s~^, compared with the errors for earlier data of ±1 km s~^. 



We collected ^-band photometry from three sources. The first was from lMarkova et al 



(l2001al ). This provided con current photo r netry for the Ha data previously published. 



The second set came from Percy et al. ( 2001 1^. These observations also ended at 



nearly the same time as the first data set. Finally, we downloaded the photoelectric 
photometry in the V^-band from the American Association of Variable Star Observers 
(AAVSO). The AAVSO data are helpful in understanding the long-term trends, and 
we only used data where measurements of the check and comparison stars differed 
from expected values by less than 0.05 mag. The errors in the AAVSO measure- 
men ts are typica l ly aro und 0.01 mag, comparable to those of Markova et al. (2001) 



and iPercy et al.l ( 120011 ). The combined set contains 3142 measurements from 1985 



to 2009. 



3. The Long-Term Photometric and Ha Equivalent Width Variability 



Markova et al.l ()2001aJ) found that the Ha line flux, obtained by correcting the 



observed equivalent widths for the changing continuum, varied in concert with the 



''Available for download at |http: / /schwab. tsuniv.edu/papers /paspc /pcyg/pcyg.html 
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V^-band flux over the period from 1989 to 1999 (see tlieir Fig. 3). Here we extend 
tlieir work by considering tfie long-term photometric and H« variations through 2007. 
Figure 1 shows the large time span of available photometry. The light curve over 
this interval shows rather modest, ±0.1 mag variations, consistent with the star's 
classification by van Genderen (2001) as a "weak-active" LBV. It exhibits the kind 
of variability associated with a short SD-phase, similar to that reported by Markova 
et al. (2001a). The short SD-phase is most evide nt in the data prior to 2000, when 



the s tar experienced two fadings of ~ 0.1 mag (IMarkova et al.l l2001at iPercy et al. 



20011 ). The amplitude of this long-term variabihty decreased in subsequent years, 
which indicates that the properties of the short SD-phase change with time. We 
made a fit of the very-long-term trend and found a brightening rate of ~ 0.17±0.01 
mag century"^ (overp lotted in Fig. 1). This rate is c onsistent with the very-long -term 



trend of 0.15 ± 0.02 mag century ^ documented by lde Groot &: LamersI (Il992[ ) 



Figure 2 presents the "dirty" discrete Fourier Transform (Roberts et al. 1987) 
of the 24.5 y l^-band photometry with the long-term brightening (Fig. 1) removed. 
There are no individual significant peaks in the periodogram, but there is a general 
tendency for more power to appear at the longer timescales (lower frequencies). Thus, 
the longer timescales of the short SD-phase variability tend to dominate the light 
curve. 

We measured Ha emission strength for both the new and originally reported 
spectra (Markova et al. 2001a) for a total of 158 measurements covering the interval 
from 1994 to 2007. Equivalent widths of the full Ha profile (including both the blue 
absorpti on and large emission component) were measured in the same manner as 



done by iMarkova et al.l ( l2001al ) in order to keep the data sets mutually consistent. 
The only improvement is that telluric H2O lines in the vicinity of Ha were removed by 
means of a template fitting procedure (telluric) in IRAF. This correction resulted 
in equivalent width increases of less than 2%. This is much smaller than the typical 
measurement error of 6% (as found by comparing equivalent widths from closely 
spaced observations. At < 2 d, where the variability of this star is minimal). Since 
the available wavelength range around H alpha does not extend beyond the electron 
scattering line wings to the actual continuum levels, a multiplicative constant was 
used to retrieve the full equivalent width of the line. This correction, WA(net) = 
1.096 ^^(Ritter), accounts for unseen line wing flux and unmeasured flux lying 
below our continuum placement (over an integration range of 6531.5 to 6593.5A) 
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and is identical to that adopted by iMarkova et al.l ( l2001al ) for the Ritter data. The 



hehocentric Juhan dates and net adjusted equivalent widths are tabulated in columns 
1 and 2 of Table 1. 

The actual line flux of Ha can be estimated by correcting the measured equiv- 
alent widths for the changing continuum. In order to make this transformation, 
we averaged the ^-band measurements made within 20 days of each spectroscopic 
measurement. This time span was chosen to cancel any small but fast variations 
and to include enough measurements for a reliable average. We compared all the 
photometry measurements to a benchmark V=4.8 to remain consistent with the flux 
correction adopted by Markova et al. (2001a). The equivalent widths were corrected 
using the relationship 



WxicoTi) = PyA(net)10-°-^(^(*)- 



-4.8)) 



The averaged V magnitudes and flux corrected equivalent widths are given in columns 
3 and 4 of Table 1. If no magnitude was available within ±20 days, then no 
correction was applied, which affects 18 of our measurements. These correction 
factors are usually small (~ 4%) and comparable to the photometric scatter within 
each time window. 

We show the temporal variations in the flux cor rected equivale n t widt hs in 



Figure 3. The plot includes earlier measurements from iMarkova et al.l ( l2001af ). the 
new measurements from Table 1, and some additional measurements from 2005 to 
2007 from Balan et al. (2010). There are two maxima (occurring around 1992 and 
2002) that are separated by ~ 10 years, which is lo nger than the reported lengths of 
the short SD-phase found by Markova et al. (j2001a ). de Groot et al. (2001), or Percy 



et al. (2001). Furthermore, the rise and fall around the peak in 1992 are steeper than 
that for the 2002 peak. There is also ample evidence of faster variability within each 
observing season that appears to be unrelated to the longer term trends. 

A visual comparison of the V^-band photometry in Figure 1 with the fiux- 
corrected Ha equivalent widths in Figure 3 immediately shows some variations in 
common. We found that the relative flux (from the time interpolated magnitude) is 
positively correlated with the corrected equivalent width. A linear fit yields a slope 
of A(F/ < F >)/A{W\/ < Wx >) = 0.16 ± 0.01, confirming the visual impression 
of co-variability. In order to compare directly the photometry and Ha equivalent 
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widths, we removed the long-term linear trend from the photometry (Fig. 1), per- 
formed a running average of the photometry differences using a Gaussian weighting 
scheme parametrized by a Gaussian FWHM, transformed the resulting flux differ- 
ences into a variation in A units according to the correlation slope given above, and 
then added the mean equivalent width to the final result. We made a number of trial 
comparisons by varying the adopted Gaussian FWHM to smooth the photometry, 
and the best fit with FWHM = 598 d is shown as a solid line in Figure 3. For com- 
pleteness, we also plot a similar running average of the Ha equivalent widths as a 
dashed line. The agreement between the temporally smoothed photometric and flux 
corrected Ha variations is striking and it appears to confirm the positive correlation 
first noted by Markova et al. (2001a). The fact that smoothing parameter values 
smaller than 598 d yield worse fits suggests that the co- variations are less correlated 
on shorter timescales. Taken at face value, this result indicates that the continuum 
and Ha emission fluxes sample structures in the wind in different ways (probably 
because of different sites of formation in the outflow). Finally, we note that the 
correlation also exists between the running averages of the continuum flux and the 
uncorrected equivalent widths W\{net), so the covariations are unrelated to the flux 
correction procedure. 



4. Ha Profile Morphology Variability 

4.1. Emission Component Changes 

The large Ha equivalent width variations described in §3 are accompanied by 
changes in the morphology of the proflle. We present two individual spectra in 
Figure 4 that represent the extrema of the equivalent widths observed (a minimum 
at HJD 2,450,004, plotted with a dashed-dot line, and a maximum at HJD 2,452,070 
plotted with a solid line). It is clear that the profile experiences a change in the 
peak emission intensity, the line width, the net profile velocity, and the shape of 
the blue absorption trough. For each of the spectra collected at Ritter Observatory 
we measured the peak intensity above continuum level, Ip, which we corrected for 
the changing continuum level in the same manner as the equivalent width (§3), the 
FWHM (profile width) of the emission portion of the profile above continuum, and 
a relative radial velocity AV^, derived by cross-correlating each profile against an 
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unweighted average of all the spectra obtained at Ritter Observatory. We chose 
to use a cross-correlation technique because this method is model-free and is most 
sensitive to the steep emission line wings, resulting in a measure similar to a FWHM 
bisector velocity. The resulting measurements of FWHM and AVr are shown for 
these two profiles in Figure 4 with horizontal and vertical lines, respectively. 

All these measurements are given in columns 5, 6, and 7 of Table 1, and are 
plotted as a function of time in the three panels of Figure 5. We see that times of 
strong emission (for example, HJD 2,452,500; see Fig. 3, MJD 5.25-10^) correspond 
to profiles with the largest peak intensity and smallest FWHM. We also see a small 
radial velocity shift that is correlated with the long-term variations. The profile had 
the largest (most positive) velocity when the line flux at the position of the emission 
peak was strongest, which was also when the profile showed the smallest FWHM 
(Fig. 5). This is hkely due to changes in the P Cygni absorption component. When 
the profile has the most emission, the blue absorption portion appears to shift to a 
more positive velocity and removes more of the blue side of the emission peak (see 
Fig. 4), and thus, the net radial velocity tends toward a larger (more positive) value 
at those times. We find that the measurement errors for A^. and the FWHM are 
about 1 km s~^, which adds in quadrature to the calibration errors discussed in §2, 
to yield net errors of approximately ±1.5 km s^^ for most of the data, and ±3.2 km 
s~^ for data taken after 2002 September. The errors for Ip are on the order of 3 — 5%. 

Kashi (2010) has suggested P Cygni is a binary system with a fainter B-type 
companion and that small long-term radial velocity variations due to reflex motion 
might be observed in extended high resolution spectroscopic observations. This can- 
not be the explanation for the AVr changes we observe, since the Ha emission is 
formed over a volume that is much larger in radius than the predicted semimajor 
axis of the putative orbit of the P Cygni primary star. Wind gas leaving the star 
at any instant would have a Keplerian orbital component that decreases with dis- 
tance from the center of mass. As the gas packet moves out to the radius where 
Ha becomes optically thin and emits the photons we observe (at ~ lOR^ and larger; 
see below), the radial outflow component will increase by radiative driving while the 
orbital motion component will drop with distance to conserve angular momentum. 
Thus, at the large distance of fine formation, the gas motion will be almost com- 
pletely radial. If the putative companion is to be found from radial velocity variations 
of this star, then detailed analyses of photospheric or wind lines formed very close 
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to the star will need to be analyzed. Further, these radial velocity variations are not 
strictly periodic, and cannot be considered orbital motion. Lastly, given the method 
of measuring these velocities, the measured radial velocity is at least partially due 
to morphological changes in the Ha line profile. 



4.2. Blue Absorption Changes 

We used all of the Ha spectra fr om Ritter Observatory, including those that were 



measured by lMarkova et al.l (j2001af ). to investigate the variations in the blue absorp- 
tion trough of the P Cygni profile. This portion of the profile is especially interesting 
as it is formed in the outflowing gas along our line of sight to the star. In order to 
emphasize the relative changes in line absorption, we first formed a reference, average 
high-intensity, minimum-absorption spectrum, as follows. At each wavelength step, 
we ordered the time-series by intensity and then constructed the average of all the in- 
tensities falling between the 90th and 95th percentiles at that wavelength step. This 
removed any spurious peaks caused by cosmic rays from contaminating the minimum 
absorption average. We then divided each of the spectra by this reference spectrum 
to form a matrix of quotient spectra. Since we are interested in the variability of 
the central absorption, and not that of the far wings, and because the line wings 
never reach the continuum in the region we recorded, the quotient spectra had a 
depressed continuum. We then re-normalized these to a unit continuum (outside of 
the velocity region ±500 km s~^). These spectra are illustrated in a gray-scale dy- 
namical spectrum in Figure 6. In this figure, we present each quotient spectrum as a 
function of radial velocity and time with a gray-scale intensity between the minimum 
value (black; 0.14 in the quotient) and maximum value (white; 1.75 in the quotient) 
based upon a linear time interpolation between the nearest observations (indicated 
by arrows). The low absorption reference spectrum is displayed for comparison in 
a panel below the gray-scale image. For simplicity, these quotient spectra were not 
corrected for the variable continuum flux since we are interested in both emission 
and absorption changes. 

We need to bear in mind that the low absorption spectrum was formed by 
different subsamples at each wavelength point, and this has important consequences 
for the appearance of the dynamical spectrum. For example, inspection of Figure 4 
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shows that blue absorption core can extend to high negative velocities (dash-dot 
line) while at other times the blue absorption is limited to moderate velocities (solid 
line). Thus, the construction of the low absorption spectrum will be dominated by 
the latter examples in the vicinity of the blue absorption edge, and in our collection 
of Ha spectra, the more extended blue absorption occurred much more frequently. 
Consequently, the quotient spectra in Figure 6 appear to be dominated by a blue 
absorption feature, near -220 km s'^ except near HJD 2,452,500 (MJD 5.25-10^) 
when the blue edge moved to a more positive velocity. This feature is due entirely 
to our selection of spectra from around HJD 2,452,500 in making the low absorption 
spectrum. 

We also see evidence in Figure 6 of several blueward moving, absorption features 
(primarily between —100 and —200 km s~^) that appear morphologically similar to 
the Discrete Absorption Components (DACs) observed previously in other spectral 
hues (Israehan & de Groot 1999; Markova 1986a, 2000). Figure 7 is a montage of a 
subset of the quotient spectra. It shows how the DAC (center) moved progressively 
blueward over this time span 800 d). There are times where the regions between 
successive DACs appear bright in the dynamical spectrum, which correspond to those 
cases (usually sparsely sampled in time) where the flux was higher than the mean in 
the 90% to 95% part of the flux distribution that defined the minimum absorption 
spectrum. Finally, we see in the gray-scale image of Figure 6 the long-term variations 
in peak intensity (Ip) and wing extension (FWHM) that are associated with the short 
SD- phase, equivalent width variations (Fig. 5). 

We measured the variability of the quotient spectra by calculating the standard 
deviation at each pixel of velocity space. This standard deviation spectrum is shown 
in Figure 8. There is a broad feature centered at rest velocity which is associated 
with the varying peak height of the profile. Another feature is present at ~ -1-150 
km s~^, which could be caused by either the variations in the profile width (FWHM) 
or red emission wing variability from traveling bumps (Markova 2000). The largest 
two features are from the DACs (seen as a broad peak around —140 km s~^) and 
from the variations present near the blue edge of the absorption core (visible as a 
peak centered at —220 km s~^). 

Figures 4, 6, 7, and 8 all demonstrate that there are significant changes observed 
in the profile near the blue edge of the absorption core. The blue-edge velocity of 
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a P Cygni profile is generally set by the position where the absorption core rises 
to intersect with the local continuum, and this velocity corresponds to the fastest 
moving gas projected against the disk of the star. However, in the case of Ha in the 
spectrum of P Cygni, this location is often poorly constrained because absorption 
may extend blueward with a shallow slope (see Fig. 4). Thus, we decided instead to 
document the kinematical changes near the blue edge by measuring the position of the 
absorption core flux minimum, ^^(min), which is normally found near —210 km s~^ 
where the slope of the profile changes sign abruptly. We determined this position 
by finding the zero crossing in the numerical derivative of a smoothed version of 
the spectrum. The S/N ratio was sufficient in all our spectra that the zero of the 
derivation was always well-defined. This estimate of the minimum flux velocity 
Vr(niin) is given in column 8 of Table 1, and the errors in V^(min) are comparable 
to the errors associated with the emission line kinematic measurements (§4.1). This 
velocity is probably related to the wind speed at a location in the wind where Ha 
ceases to be optically thick. 

It is difficult to measure the radial velocities of the DACs because their mor- 
phologies vary and because the absorption may consist of multiple components. We 
decided to measure a centroid for the DACs wherever possible by means of the rela- 
tionship 

vJl — Q(Vr)) dVr 
jy, (1 - Q{Vr)) dVr 

where Q represents the quotient spectrum and Vr is the radial velocity. We adopted 
a velocity range of Vi — —200 km s~^ and V2 — —100 km s~^ based upon the 
strongest regions of the standard deviation spectrum plotted in Figure 8. Typical 
errors in these measurements (from the scatter in densely sampled regions of the 
time series) are ±3 km s"^. This approach worked for most of the spectra, but some 
low contrast features were not measured correctly, and are omitted from Table 1 and 
our analysis. During some epochs, there were multiple DACs present, so K(DAC) 
represents a weighted average of multiple components. The DAC radial velocities 
Vr{DAC) are given in column 9 of Table 1. Column 10 lists a relative equivalent 
width for the DAC measured by a direct numerical integration of Q{vr) between 
vi and V2- The typical errors associated with these equivalent width measurements 
are on the order of 5%, which is similar to the errors associated with the equivalent 
widths of the profile (§3). 
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We show the temporal evolution of the DAC and blue-minimum flux velocities 
in Figure 9. At some epochs, a long progression of DAC velocities is present. The 
most well defined sequence began near HJD 2,450,700, and lasted about 800 days 
(Fig. 7). This timescale is much longer than that of a typical DAC observed in 
an O star, where the progression is measured in hours or days. As this particular 
DAC was relatively narrow and recorded in many spectra over its duration, it is an 
optimal feature to measure the Ha DAC acceleration. From a simple linear fit, shown 
overplotted in Figures 7 and 9, we measured the acceleration to be —0.047 ± 0.002 
km s~"^ d~^. For comparison, the observed DACs in the spectrum of a normal hot 
supergiant of similar spectral type, e Ori (HD 37128; BO la), have an acceleration of 
-500 km s-^ d-^ (Prinja et al. 2002). 

We used the Lomb-Scargle periodogram method (Scargle 1982) to search for 
a recurrence time in the appearance of the DACs, and we derived a cycle time of 
1700 d. This recurrence timescale is shown in Figure 9 where we overplot the hnear 
acceleration of the DAC shown in Figure 7 for three additional epochs over the time 
span of the Ritter data. There is some evidence that a DAC progression is seen at 
each of these four epochs. The major deviations from the expected velocity trends 
occur when multiple components are present. For example, there were two DACs 
present between HJD 2,451,700 and 2,452,400, and the velocity centroid we measured 
represents a blend of these components. Neither of these two DACs occurred at the 
expected recurrence time in the 1700 d cycle. 

Our work represents the first detection of DACs in the blue absorption trough of 
Ha, and their properties differ from those observed in other spectral lines (Israelian 
et al. 1996; Markova 1986a, 2000). For example, the recurrence timescale of 1700 d 
is much larger than the 200 d interval found in earlier work, and the acceleration we 
measure is about a factor of 10 smaller than measured by others. We suspect that 
these differences are due to the large optical depth of the Ha line compared to that 
of other lines where DACs have been investigated. This will mean that the radius 
of optical depth unity is larger for Ha (Najarro et al. 1997), and consequently any 
wind structures formed at smaller radii will have no affect on the Ha line formation. 
We suspect that observations of other, less optically thick lines are more sensitive to 
the detection of DACs formed at smaller radii in the wind of P Cygni where more 
and faster accelerating structures may exist. 
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5. Discussion 

The K-band and H« variations we observe need to be interpreted in the context 
of current models for the star and its wind. Langer et al. (1994) discuss the atmo- 
spheric properties of P Cyg, and they argue that the star is in the LBV phase where 
the temperature and hehum abundance are increasing, and the mass and luminosity 
are decreasing, as the star evolves towards the Wolf-Rayet phase. Langer et al. em- 
phasize the earlier conclusion from Pauldrach & Puis (1990) that the atmospheric 
parameters are close to a bi-stability point where the mass loss rate can change by 
an order of magnitude with small changes in radius and/or luminosity, which may 
explain the great eruptions observed in prior centuries. The atmospheric parameters 
are well established through a detailed quantitative spectroscopic analysis by Najarro 
et al. (1997) and Najarro (2001), who find that He is overabundant and that the mass 
loss rate is high (f» 2 x 10^^ Mq y~^ including wind clumping effects) and wind ter- 
minal velocity is low [v^o ~ 185 km s"-*^). Najarro et al. (1997) derive a systemic 
velocity of 7 = —29 km s~^, and thus our minimum measurement of V^in = —215 
km s~^ is consistent with their estimate of 7 — i^oo = —214 km s~^ as this velocity 
measurement is related to Voo- They estimate that the continuum forming radius is 
76Rq, which for a distance of 1.8 kpc implies an angular diameter of ^ = 0.39 mas. 
On the other hand, Najarro et al. (1997) predict that the emitting size of H« will 
be much larger because of its greater optical depth. For example, their models show 
that there is a local maximum in the wind temperature distribution (presumably 
where the recombination processes that form Ha peak) near r/R^ — 11 (see their 
Fig. 5b). The corresponding angular size for Ha of ~ 4 mas agrees well with the 
range of 3 — 7 mas from Ha interferometry by Balan et al. (2010). Thus, we need to 
keep in mind that the Ha variations reflect changes over a much larger spatial scale 
in the wind than those observed in the y-band flux. 

Variations in the Ha emission equivalent width are related to changes in both 
the mass loss rate and the wind velocity. In a very simplifled approach, we can 
assume that most of the Ha flux originates in the optically thick region projected on 
the sky, 

/ = 2nrlF{T) 

were r-^ is the boundary separating the optically thick and thin regimes, F{T) is 
the monochromatic surface flux, and T is the wind temperature at rr (Najarro et 
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al. 1997). If we assume that the wind is approximately isothermal at this physical 
location (a reasonable choice: see Fig. 5b in Najarro et al. 1997), then the emission 
flux variations are due to changes in the projected size of the optically thick region, 

A/// = 2Ar,/r,. 

Thus, we expect that the relative variations in angular size will be only half as large 
as the emission equivalent width variations, which is probably consistent with the 
lack of measurable size changes in the Ha interferometric measurements (Balan et 
al. 2010). 

The Ha optical depth is dependent on the electron density squared since the 
emission is a recombination process. Thus, we expect that the optical depth unity 
boundary will always be defined by the location in the wind with a specific char- 
acteristic density, p^-- We assume that Pr has an approximately constant value so 
that the effective boundary will vary as fluctuations in the wind mass loss rate 
and velocity define the radius where the density reaches pr- According to the mass 
continuity equation, Vr is related to this density by 

„ M 

r = 

ifrprV 

where M is the mass loss rate and v is the wind velocity at the radial distance r^- 
We can differentiate the mass continuity equation to express the radius variation in 
terms of the changes in M and v, 

2rrArr = -^A[M/v], 
inpr 

which we divide by to obtain 

^Arr ^ A[M/v] 

Since we argued above that the flux also varies as r^, we can then use the relation 
above to re- write the fractional flux variation in terms of logarithmic changes in M 
and 

Aln/ = AlnM- Aln^;. 
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Since we have observational data on the variations in emission strength and wind 
velocity, we can rearrange this equation to solve for the mass loss variations as a 
function of flux variations, 

AlnM/Aln/ = 1 + Alnv/Aln/. 

Pills ct al. (1996) present a much more detailed analysis of the dependence of the 
emission equivalent width on the wind parameters of hot, massive stars. However, 
in the limit of high optical depth, their expression for the emission flux (their eq. 41) 
leads to a similar relation, 

AlnM/Aln/ = ^(1 + Alni;/Aln/). 

We found in the previous section that the Ha equivalent width appears to vary 
inversely with two quantities related to wind dynamics, the Ha emission peak FWHM 
and the blue minimum flux velocity V^(min). Figure 10 quantifies this relation- 
ship. The upper panel shows the inverse correlation between the emission peak 
FWHM and fiux corrected Ha equivalent width, and if we take FWHM as a proxy 
for the wind speed, then a linear fit of natural logarithms of these measures gives 
A In t;/ A In/ = -0.66 ± 0.12. We caution that the FWHM is also influenced by 
the absorption component of Ha, and we showed above (Fig. 4) that the absorp- 
tion component moves inward towards the line core when the emission is strong. 
Consequently, the apparent decrease in FWHM as the emission increases probably 
results both from a wind speed decrease and a blue wing decline due to blending 
with the absorption component. The minimum flux velocity is perhaps a more di- 
rect measurement of wind speed (at least in our line of sight), and we show in the 
lower panel of Figure 10 the co- variations of the difference between V^(min) and the 
systemic velocity of P Cyg, 7 = —29 km s^-*^ (Najarro et al. 1997), as a function of 
the corrected equivalent width. A fit of the logarithmic slope here gives a smaller 
estimate of A In t;/ A In / = -0.22 ± 0.04. 

If we adopt the minimum flux co-variation result as representative of the wind 
velocity component of variability, then the mass loss rate variation we derive from the 
relation above has an emission flux dependence of A In M 0.78 A In /. Omitting the 
bottom and top 10% of the distribution of Wx{coyy), the derived range in emission 
strength in our observations of ±14% probably implies mass loss rate changes of 
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±11%. Markova et al. (2001a) used the optically thick relation from Puis et al. 
(1996) to arrive at an estimate of ±9% for the mass-loss variation amphtude. We 
showed above that the relation from Puis et al. carries a factor of 3/4 that is missing 
from our simple analysis, and if we use the Puis et al. relation instead, we arrive 
at a mass loss rate variation of ±8%, confirming the estimate from Markova et al. 
(2001a). In this scenario, the Ha emission variations result from changes in the 
effective emission radius caused by variations in the mass loss rate and wind 
velocity. During episodes when the mass loss rate is higher and the wind is slower, 
the projected size of the emission region increases leading to larger Ha emission flux. 
The same process probably causes the y-band variations, but the fractional radius 
variations must be smaller at the continuum forming radius because we found in last 
section that Aln/[F]/Aln/[HQ;] = 0.16 so that Alnr^-y = 0.16A InrT-Ha, i.e., the 
continuum size variations are only 16% as large as those in Ha. It is possible that 
the changes in the mass-loss rate are caused by a change in the luminosity of the 
star which would propagate through the wind, and be observed in both the V-band 
brightness and the emission line flux of Ha. 

Finally we return to the relationship of the DACs to the short SD-phase vari- 
ations. We found a trend with the DAC strength and the Ha emission equivalent 
width. We show in Figure 11 the temporal behavior of DAC quotient equivalent 
width 1^;^(DAC) (Table 1, column 10) along with scaled, running averages of the Ha 
equivalent width and F-band flux (Fig. 3). We rescaled the amphtude of the Ha flux 
by AWx(DAC) — Al^A(corr)/20.46 and the photometric hght curve was rescaled by 
AiyA(DAC) = Afv X 22.56. These curves were then shifted vertically to match the 
Wa(DAC) points. We see that the DAC strength variations track both the Ha and 
y-band fiux variations, suggesting that the DACs are related in some way to the 
short SD-phase changes. For example, we see that some of the strongest DACs were 
observed when the Ha emission was strong (around HJD 2,452,500; MJD 5.25-10^ 
in our plots) and the DACs were very faint or absent when the emission was weak 
(around HJD 2,450,500). 

The DAC phenomenon is primarily observed in the UV wind lines of hot stars 
(Kaper et al. 1999; Puis et al. 2008), and, in fact, DACs are only rarely seen in the 
Ha profile where wind-related variations are usually due to changes in the dense and 
slower moving wind close to the star (Kaper et al. 1997; Markova et al. 2005). The 
only comparable DAC observed in Ha was in HD 92207 (Kaufer et al. 1996) with a 
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lifetime of ~ 150 d and only one DAC observed, so the recurrence time is unknown. 
This star is cooler than P Cygni and is a "normal" supergiant, compared to the 
luminous blue variable nature of P Cyg. The DACs observed in the UV wind lines of 
0-stars are first seen at velocities of 0.2 to 0.4 v^o and then they migrate blueward to 
on time scales of a day or so, exhibiting accelerations that are much slower than 
expected for the wind velocity law (Kapcr ct al. 1999). Many of these same features 
are seen in the DACs in Ha for P Cygni, although they occur on vastly longer time 
scales. For example, the wind flow time scales as i?*/foo, and while the wind gas will 
accelerate from O.lvoo to O.Qvoo in 0.5 d for an 0-star like ^ Per (Kaper et al. 1999), 
it will take some 43 d for the wind of P Cygni. However, this flow time is very short 
compared to the longevity of the DACs (10^ d), so we are led to the same conclusion 
found for the 0-stars, namely that the DACs represent some kind of perturbation in 
the wind through which the gas flows. Changes in wind velocity and/or mass loss 
rate can cause shocks and create structure in the wind, and these structures produce 
density enhancements and/or velocity plateaus that imprint DACs in the wind lines 
(FuUerton & Owocki 1992; Runacres & Owocki 2002; Puis et al. 2008). 

Current theory suggests that the DACs are related to large scale spiral features in 
the wind known as co-rotating interaction regions that are formed at the intersection 
of fast and slower outflows, which develop from some inhomogeneity in the stellar 
photosphere (for example, pulsation or spots; Cranmer & Owocki 1996; Lobel & 
Blomme 2008). In these models, it is the slow transit of these equatorially centered 
regions across the photosphere that creates the DACs in the absorption cores but has 
little influence on the emission parts of the wind fine (Dessart 2004; Lobel & Blomme 
2008). However, in the case of P Cygni, we flnd that emission parts do appear to 
strengthen when DACs are prominent (Fig. 11), and this indicates that the wind 
perturbation profoundly affects both wind gas surrounding the star and the wind gas 
projected against the photosphere. Thus, we suggest that the structures causing the 
DACs in P Cygni may be more spherically symmetric than assumed in the geometry 
of the co-rotating interaction regions. In some models the seed perturbation occurs 
at a fixed longitude on the star, so that a new wind structure appears each time the 
star rotates (although for the best studied case of HD 64760, Lobel & Blomme 2008 
argue that the originating spots must rotate some five times slower than the star in 
order to fit models to the observations). If this is the case for P Cygni, then the 
1700 d DAC recurrence time may be related to the star's rotational period. Markova 
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(2000) found an upper limit of 100 d for the rotational period based upon estimates of 
the stellar radius and the projected rotational velocity. However, the line broadening 

in early supergiants may be dominated by turbulence rather than rotation (Howarth 
et al. 1997, 2007; Markova & Puis 2008), so a longer rotational period remains a 
possibility. However, regardless of the origin of the DACs, their close relation with 
the emission line and continuum flux variations (Fig. 11) suggests that much of the 
short-SD variability is caused by propagating structural perturbations in the outer 
atmosphere of the star. 

The discovery of the short SD-phasc variations in P Cygni and its relationship 
to the wind velocity and DAC occurrence in the wind is a new observational result 
that warrants future investigation both for this star and other LBVs. The changing 
characteristics over these long timescales may eventually lead to a better under- 
standing of the LBV stage of evolution and the underlying physics of their winds 
and circumstellar environments. We are currently pursuing a three year, spectro- 
scopic monitoring program of Galactic and Magellanic Cloud LBVs. Such long-term 
observations will reveal if DACs are present in all LBVs and will show whether the 
variations found in P Cygni are a general phenomena among LBVs. 
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Fig. 1. — The V-band variability of P Cygni between 1985 and 2009. The long-term 
changes are representative of the short SD-phase. We overplotted the fit for the 
very-long-term brightening of the star with a solid line. 
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Fig. 2. — Discrete Fourier transforms of the l/-band photometry. The bottom panel 
shows a close up of the long timescale region that is associated with the short-SD 
phase and 100 d timescales of variabihty. The overplotted stair steps show the 
amplitudes rebinned into increments of 0.001 cycles d~^. 



-35 - 




Fig. 3. — A direct comparison of the smoothed F-band photometry (sohd hne) and 
the flux corrected Ha equivalent widths. The photometric hght curve is a running 
average of the F-band flux that was re-scaled and offset to match the equivalent 
width variations (see text), with the scale of the differential light curve given on 
the secondary y-axis. A running average of the corrected equivalent width mea- 
surements is also overplotted as a dashed line. The equivalent width measurements 
from Markova ct al. (2001a) arc denoted by diamonds, our new measurements arc 
represented as triangles, and the measurements of Balan et al. (2010) by squares. 
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Fig. 4. — A comparison of two observed Ha profiles corresponding to extremes of 
the emission variability (maximum and minimum in corrected equivalent width). 
The spectrum plotted as a solid line is from HJD 2,452,070, while that shown as a 
dash-dot line is from HJD 2,450,004. Vertical and horizontal hne segments show the 
velocity offset AVr and the FWHM range, respectively. 
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Fig. 5. — The long-term variability of the FWHM (top panel), the relative radial 
velocity AV^ (middle panel), and flux corrected peak intensity Ip (bottom panel) of 
the emission peak of the Ha profile. The long timescales of variability are similar to 
those seen in the F-band photometry and the Ha equivalent widths (Fig. 3). 
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Fig. 6. — A dynamical spectrum showing the quotient of the Ha spectra and a 
constructed low absorption profile (shown in the lower panel) as a function of radial 
velocity and time. The blueward moving features are the DACs, which are almost 
always present in Ha. The gray-scale image is formed by a time interpolation (over 
100 d in some rare cases) to fill in observational gaps in the data. Arrows on the 
right of the diagram indicate the dates of observation. The range of the plot is from 
0.14 (black) to 1.75 (white) in intensity ratio. The original spectra, as well as the 
maximum intensity profile, are available in a supplementary file in the online version. 



-39- 




Fig. 7. — A series of quotient spectra offset according to the time of observation 
(HJD — 2,400,000 indicated for the first and last spectrum in this sequence). The 
DAC present migrated from approximately —125 km s^^ to —160 km s^^ over this 
interval. The dotted line represents a linear fit to the measured centroid velocities 
for this DAC. 




Fig. 8. — The pixel- by-pixel standard deviation of the quotient spectra (sohd line). 
We also show the low absorption reference profile (dotted line; re-scaled to this 
range) to highhght those parts of the profile where the largest relative variations are 
occurring: near the terminal velocity blue edge (near —220 km s ~^), over the range 
traversed by the DACs (centered near —150 km s ~^), near the emission peak (0 
km s ~^), and in the emission wings (±150 km s ~^). The standard deviation of the 
quotient is larger in the absorption core because the low absorption profile is close 
to zero there. 
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Fig. 9. — The temporal variations of the measured velocities of the DACs, Vr{T)AC) 
(plus signs), and the minimum flux wind velocities, K-(niin) (diamonds). The dashed 
hne represents the acceleration fit to the DAC progression near HJD 2,451,000, and 
the dotted fines represent the same fit translated by intervals of 1700 days, which we 
derived as a possible recurrence time. 
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Fig. 10. — A comparison between the flux corrected equivalent width of Ha and 
the FWHM of the emission peak (top panel). The best fit (solid line) has a slope 
of —0.66. The lower panel shows a comparison between the corrected equivalent 
width of Ha and the difference between the minimum flux velocity V^(min) and the 
systemic velocity 7 of P Cyg. A linear fit yields a slope of —0.22 (solid hue). 
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Fig. 11. — The temporal variations of the relative equivalent width of the DACs. 
We overplot the running average of the V-band photometry (solid line) and the Ha 
corrected equivalent width (dashed hne) that were both rescaled and shifted to match 
the DAC variability (see text). The DAC strength appears to follow the variations 
associated with the short-SD phase. 



